Rheology and Contact Lifetime Distribution in Dense Granular Flows 
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We study the rheology and distribution of interparticle contact lifetimes for gravity-driven, dense 
granular flows of non-cohesive particles down an inclined plane using large-scale, three dimensional, 
granular dynamics simulations. Rather than observing a large number of long-lived contacts as 
might be expected for dense flows, brief binary collisions predominate. In the hard particle limit, 
the rheology conforms to Bagnold scaling, where the shear stress is quadratic in the strain rate. As 
the particles are made softer, however, we find significant deviations from Bagnold rheology; the 
material flows more like a viscous fluid. We attribute this change in the collective rheology of the 
material to subtle changes in the contact lifetime distribution involving the increasing lifetime and 
number of the long-lived contacts in the softer particle systems. 
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The rheology of granular materials is relevant to many 
areas of nature and industry, from mountain avalanches 
and mud-slides, to grain transport and storage [l], 0- 
We study a particularly simple type of flow reminis- 
cent of avalanches: gravity-driven, dense granular flow 
down a rough, inclined plane. This geometry is the 
archetypal granular flow with which one can study the 
relation between the stress state and the dynamics and 
structure, i.e. the constitutive relation of the granular 
material. Indeed, a number of recent, well-controlled 
experiments [H, [j, @, and large-scale simulations 
0, 0, B E3] have motivated numerous theories that 
capture some of the features of inclined plane flows 
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To date, however, the most generally accepted treat- 
ment of granular rheology is still the physical picture put 
forward by Bagnold over 50 years ago [20|,[21|. Bagnold 
described a mechanism of momentum transfer between 
particles in adjacent layers that assumes instantaneous 
binary collisions between the particles during the flow. 
Under this assumption, the inverse strain rate is the only 
relevant time scale in the problem leading to a consti- 
tutive relationship between the shear stress a and strain 
rate 7 of the form 



(1) 



where k is independent of 7. 

Despite recent concerns regarding the validity of Bag- 
nold's original experiments [22|, and the applicability of 
the theory [23j], Eq. Q] has proved rather successful in 
predicting the rheology of dense granular flow down an 
inclined plane 0, Q. This is somewhat at odds with 
the original understanding of the Bagnold theory as it 
is based on the physics of rapidly flowing hard-spheres 
where binary collisions predominate. Dense flows [2J], on 
the other hand, are thought to be controlled by enduring, 
multiple interparticle contacts forming extended stress- 
bearing structures [lH and/or large length-scale cooper- 



ative dynamics In fact, the stresses arising 

from such contact forces in dense flows are practically 
an order of magnitude larger than the kinetic stresses 
associated with velocity fluctuations |8|. 

In this letter, we analyse both the constitutive relation 
and the interparticle contact dynamics of dense flows of 
non-cohesive granular particles. To discuss interparticle 
dynamics we introduce two timescales. The first is the 
lifetime of interparticle contacts and the second is the 
inverse strain rate in the material. The inverse strain rate 
or shear time sets the fundamental timescale over which 
particle rearrangement events occur during the flow. 

We find that the typical contact lifetimes in the system 
are of the same order of magnitude as the binary colli- 
sion time and insensitive to the location of the pair of 
contacting particles within the flowing pile. In contrast, 
the shear time is strongly dependent on the height and 
varies considerably over the parameter space of interpar- 
ticle interactions studied here. The scenario that emerges 
is somewhat counter-intuitive: Even in dense flows the 
dynamics at the microstructural scale remain predomi- 
nately that of short-time, binary collisions. When the 
interparticle contact lifetime is significantly shorter than 
the shear time, the system is effectively in the hard- 
particle limit, and the Bagnold constitutive relation, Eq. 
1, holds. Remarkably, the Bagnold relation is extremely 
sensitive to the development of a small population of 
long-lived contacts on the timescale of the shear time. 
Changes in the interparticle interaction that lead to the 
enhancement of these longer-lived contacts also lead to 
the breakdown of the Bagnold relation. In such cases it 
must be supplemented by an additional term linear in the 
strain rate. Thus the functional form of the constitutive 
law is sensitively dependent on changes in the long-time 
tail of the interparticle contact time distribution. 

Our simulations are based on the model developed 
by Cundall and Strack j^, and Walton [Hj], and has 
been shown to quantitatively match experimental data 



BBS- We study N = 35, 900 monodisperse spheres 
of diameter d and mass m, flowing on a rough base of 
length 20d, width 20c?, and tilted an angle 9 with respect 
to gravity. We use periodic boundary conditions in the 
flow and vorticity directions. The height h, of the flowing 
pile is between 90d < h < lOOd depending on the angle 
of inclination. For most of the results presented here, the 
inelastic collisions are modelled as a Hookean spring and 
dashpot for forces both normal (n) and tangential (t) to 
the interparticle contact plane, but similar results were 
found for Hertzian contacts [30( . At existing contacts we 
include a friction model that obeys the Coulomb yield 
criterion with a coefficient of friction [i = 0.5. 

For Hookean contacts the coefficient of restitution pa- 
rameterises the dissipative nature of the interparticle col- 
lisions and is given by e n = exp (— 7 n i C oi/2). Here t co \ is 
the binary collision time, 

i col = ^[2fc n /m- 7 2/4]- 1/2 , (2) 

where fc n and 7„ parameterise the elastic and dissipative 
normal interparticle forces respectively. By simultane- 
ously adjusting these parameters, we fix e n = 0.88 for 
the results presented here. Varying e n has little effect 
on the lifetime distributions and the conclusions of this 
work. For the tangential interactions we set fc t — 2k n /7 
and 7t = 0. 

We explore a range of spring constants fc n , from 10 3 - 
10 9 mg/d. For fc n < 2 x 10 5 mg/d, the time step for nu- 
merical integrations St ~ 10~ 4 t, where time is in units 
of t = y/d/g. St is adjusted accordingly for larger fc n , 
to ensure accurate dynamics. After reaching a dynam- 
ical steady state, characterised by time-independent to- 
tal energy, where the mean flow velocity ranged from 
10-lQOs/dg depending on the system parameters, we col- 
lected data over a period of ~ 10 6 St. 

In Fig.QJa) we show the distributions, P(t* = T c /t co \), 
of two-particle contact lifetimes, r c , normalised by the 
binary collision time t co \, for different k n . Under this 
time rescaling all distributions exhibit a prominent short- 
time peak near the binary collision time r* — 1, and an 
approximately exponential decay towards longer contact 
lifetimes. As k n is reduced this exponential tail becomes 
broader indicating an increasing density of enduring con- 
tacts. The data in Fig. [TJis depth-averaged over contacts 
in the flowing pile away from the bottom boundary and 
the top saltating layer, at an inclination angle, 6 = 23°. 
The distributions are qualitatively unchanged over the 
range of accessible angles. We return to the more subtle 
dependence of r c on 9 in our discussion of Fig. [5] 

The normalised mean contact lifetime exhibits essen- 
tially no depth dependence as shown in Fig. H^a) . Given 
the stress-free boundary condition at the free surface of 
the pile, 7 — > there so the normalised inverse strain 
rate, 7* = 7 _1 /*col) shown in Fig.[2^b), must be height 
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FIG. 1: Dependence of the distributions P(t*) of the depth- 
averaged normalised contact lifetimes, t£ = T c /t co \, on parti- 
cle stiffness k n = 2 x [10 9 , 10 7 , 10 5 , 10 4 , W 3 ]mg/d. Data for 
9 = 23°, but the distributions are only weakly sensitive to 6. 





dependent. It also depends strongly on k n 
several orders of magnitude in our data set. 



varying by 



FIG. 2: Depth profiles of, (a) the normalised contact time 
t* (line is a guide to the eye) and, (b) the normalised shear 
time 7* = 7 _1 /t C oi, for different particle stiffness, k n — 2 x 
[10 9 (A), 10 7 (V), 10 6 (>), 10 5 (0), 10 4 (□), 10 3 (o)]m 5 /d. 
z = defines the bottom of the pile. Flow at = 23°. 



As shown in Fig. [3l over several orders of magni- 
tude in fc n , the mean normalised contact time remains 
nearly constant while the maximum contact time T Cmax , 
extracted from the distributions in FigJ.lj decreases with 
increasing stiffness, reflecting the narrowing of the con- 
tact time distributions as the particles become harder. 
Additionally, the average number of contacting neigh- 
bours per particle n c , shown in Fig.[3Jb), tends to the bi- 
nary collision limit (= 0) as the particles become harder. 
Thus, from a coordination point of view the flows are 
not densely packed. Taken in combination, the data pre- 
sented in Figs. [T]-[3]show that increasing the stiffness of 
the interparticle interaction reduces the weight of the ex- 
ponentially small long-time tail of the contact lifetime 
distribution. This narrowing of the contact lifetime dis- 
tribution is not observable from the more coarse measure 
of the mean contact lifetime distribution, which remains 
close to the binary collision time i co i. Thus the parti- 
cle dynamics are always dominated by binary collisions, 
but changing the effective particle stiffness causes subtle 
changes in the number of rare, enduring contacts. 
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FIG. 3: Dependence on the particle stiffness k n , of, (a) the 
depth-averaged contact lifetime < r * > (filled circles) and 
maximum contact time r* max (open squares), both normalised 
by tool i and (b) the depth-averaged number of contact neigh- 
bours per particle, n c . Flow at = 23°. 
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Since the interparticle contact lifetime is gcncrically 
small compared to the shear time, one might imagine 
that momentum transport is dominated by effectively in- 
stantaneous binary collisions, and might expect that the 
Bagnold constitutive law holds. To test this we now turn 
to a characterisation of the granular rheology by fitting 
the velocity profile of the flowing material to the predic- 
tion of a modified Bagnold relation of the form 



a = K7 + f3j, 



(3) 



where the coefficients k and (3 are determined by a least- 
squares fit to the velocity data. To characterise the de- 
viations from standard Bagnold rheology we define the 
ratio of the shear stress in the Bagnold and linear forms 
to be ft = ft/nj. 

Figure 2] shows the dependence of fl on k n for 9 = 23°. 
In the hard-particle limit we expect £1 — > 0. This is 
practically achieved already for fc n > 2 x 10 5 mg/d, indi- 
cating that this value of k n is appropriate for modelling 
hard (glass) particles in a system of this size 0, E|. As 
the particles are made softer, il grows meaning that the 
constitutive law approaches a linear relation reminiscent 
of a viscous fluid. Since SI is inversely dependent on j, 
it is no surprise that it grows monotonically as one ap- 
proaches the free surface as is shown by the inset to Fig.|U 
We have examined the (weaker) dependence of f2 on all 
other particle interaction parameters, e.g. friction and 
inelasticity; these results will be discussed elsewhere [3(| • 

For these dense flows the contact dynamics of all sys- 
tems studied are dominated by binary collisions, that en- 
dure for no more than t co i, such that (r c ) A f <C 1. But their 
rheology changes dramatically with fc n , going from the 
quadratic-in-shear-rate Bagnold law for hard particles to 
a nearly linear or viscous relation for softer particles. The 
dramatic change in the constitutive relation is controlled 
by a more subtle feature of the contact time distribution 
than simply the mean value. The significant rheological 
change appears to be due to the growth of the lifetime 
and number of long-lived contacts in the softer systems. 



FIG. 4: Dependence of the parameter Q on the stiffness k n 
for = 23°. Taken from the middle of the pile at h/2. Inset: 
Height dependence of f2, for k n = 2 x 10 4 mg/d. 



We emphasise that the fraction of such enduring con- 
tacts remains small in comparison to the more common 
short-lived contacts with r c ~ t co \. 

The growth of a small population of long-lived contacts 
leads to the formation of transient stress-bearing struc- 
tures within the flowing material, which can be Theolog- 
ically significant. These structures span streamlines in 
the flow and thus elastically transmit stress across their 
length in proportion to the rate of particle impacts with 
these structures ~ 7. This picture is also consistent with 
the rapid-quasi-static transition in shear flows 3l| . This 
reasoning suggests that the shear rate is the appropri- 
ate local clock with which to measure contact lifetimes 
t c so that the dimensionless contact lifetime r c 7 is the 
fundamental quantity controlling deviations from Bag- 
nold rheology in the flowing state. Given similar contact 
time distributions, faster flows should then deviate more 
strongly from the Bagnold prediction than slower flows. 

To explore this point, we use the inclination angle of 
the pile to adjust the overall shear rate and study the re- 
sulting change in the observed rheology as parameterised 
by il. In Fig. [5] we plot vs 717 for soft flows with pa- 
rameters spanning our range of stable angles and particle 
stiffnesses. Here, 77 is the lifetime of long-lived contacts in 
the system, as defined by the criterion that these contacts 
are sufficiently far out in the long-time tail of the nor- 
malised contact-time distribution so that P{n) = 10~ 3 . 
For these systems, where we see pronounced deviations 
from the Bagnold constitutive relation, £1 increases with 
increasing angle as shown by the inset to Fig.[5j All vari- 
ations of f2 can be collapsed onto a single master curve 
showing that the system's rheology depends (exponen- 
tially) on the dimensionless contact lifetime 777. The 
increasing value of these longer-lived contacts measured 
relative to the fundamental clock-rate 7 drives the sys- 
tem from Bagnold to viscous rheology. Simulations with 
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f2 consistent with or oo are omitted from Fig[5l but fol- 
low the displayed trend. Attempting to similarly collapse 
all data sets using the dimensionless mean contact time 
j(r c ) fails. Only the change in the long-time tail of the 
contact distribution correlates with the observed changes 
in the rheology. 
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FIG. 5: Plot of Q vs. 777 as defined in the text for Hookean (o: 
k n = {10 3 , W 4 }mg/d) and Hertzian (o: fc n = {10 4 , 10 5 }mg/d) 
systems at inclination angles shown in the inset, which shows 
the dependence of Q on 8 on a linear scale. The inclination 
angle is the most sensitive control of the shear rate 7. 

In summary, we have provided extensive simulation re- 
sults on the rheology and interparticle contact statistics 
of gravity-driven, dense, granular flows of non-cohesive 
grains. We observe a transition from a Bagnold consti- 
tutive relation to one reminiscent of a Newtonian fluid 
as the particles are made softer. Based on our numeri- 
cal data, and to account for this transition, we suggest a 
generalised Bagnold relation [5§| to quantify the chang- 
ing rheology. Furthermore, despite the naive guess that 
the flow in such dense systems should predominantly in- 
volve long-lived, stress-bearing structures, it turns out 
that for both hard and soft particles the microscopic par- 
ticle dynamics is dominated by frequent, short-time, bi- 
nary collisions. When examining the entire contact life- 
time distribution, however, we note that as the particles 
are made softer there is a growth in the width of the 
distribution and that the lifetime of rare, long-lived con- 
tacts grows relative to the shear time. We propose that 
the emergence of these atypically long-lived contacts is 
related to the dramatic change in the granular rheology 
with particle stiffness and thereby rationalise the surpris- 
ing result that larger inclination angles, and hence faster 
flows are actually less Bagnold in rheology than slower 
flows at smaller inclination angles. It appears that the 
constitutive relation of granular media interpolates be- 
tween these extremes and is controlled by a combination 
of the particle hardness and flow rate. 
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